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ABSTRACT 

We carry out numerical simulations of dissipationless major mergers of elliptical galax- 
ies using initial galaxy models that consist of a dark matter halo and a stellar bulge 
with properties consistent with the observed fundamental plane. By varying the den- 
sity profile of the dark matter halo (standard NFW versus adiabatically contracted 
NFW), the global stellar to dark matter mass ratio, and the orbit of the merging 
galaxies, we are able to assess the impact of each of these factors on the structure 
of the merger remnant. Our results indicate that the properties of the remnant bulge 
depend primarily on the angular momentum and energy of the orbit; for a cosmo- 
logically motivated orbit, the effective radius and velocity dispersion of the remnant 
bulge remain approximately on the fundamental plane. This indicates that the ob- 
served properties of elliptical galaxies are consistent with significant growth via late 
dissipationless mergers. We also find that the dark matter fraction within the effective 
radius of our remnants increases after the merger, consistent with the hypothesis that 
the tilt of the fundamental plane from the virial theorem is due to a varying dark 
matter fraction as a function of galaxy mass. 

Key words: galaxies: fundamental parameters - galaxies: structure - galaxies: evo- 
lution - dark matter - methods: TV-body simulations 



1 INTRODUCTION 

Galaxy merging is thought to be the dominant pro- 
cess in the formation and evolu tion of elliptical galax- 
ies (e.g.. iTbomre fc Toomrel ll972'l, This picture is broadly 
supported by both observations and simulations and fits 
naturally into the ACDM hierarchical cosmology (e.g., 
ISteinmetz fc Navarro! 120021) . Much, however, remains to 
be understood about the gas and stellar dynamics dur- 
ing mergers. For example, recent simulations indicate that 
not all mergers of ga s-rich spiral galaxies result in a bulge- 
dominated remnant (ISpringel fc Hernauisdliool . 

Although gas-rich starbursts may form many of the 
stars contained in elliptical galaxies, it is also likely that 
stars are later assembled into larger remnants via nearly 
dissipationless merg ers, particularly during the formation of 
groups and clusters \Gao et al.ll2004lh There is in fact direct 
observational evidence for "red" mergers (lacking significant 
star formation), both in local observations with the Sloan 
Digital Sky Survey (SDSS; Masjedi and Hogg, private com- 
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munication) and in a lum inous X-ray cluster at z = 0.83 
llvan Dokkum_et_aL| 1999). In addit ion, the semi-analytic 
study of iKhochfar fc Bu rkert (2003) found that the most 
recent merger of a luminous elliptical galaxy {Mb < —21) 
is likely to have been between two bulge-dominated galax- 
ies rather than two gas-rich spirals, again pointing to the 
importance of dissipationless merging for the evolution of 
elliptical galaxies. 

Observationally, early type galaxies exhibit a well- 
defined correlation among their effective radii R e , luminosi- 
ties L (or equivalently surface brightness / oc L/R 2 ), and 
central velocity dispersions a. This "fundam ental plane" 
iPiorgovski fc Davislll987t iDressler et alJll98^) is often ex- 
pressed as 

R e oc a a I b (1) 

and has been meas ured by many groups with vary i ng re- 
sults. For example, Ijorgensen. Franx. fc Kiaergaardl Jl996T) 
measured a = 1.24 ± 0.07 and b = -0.82 ± 0.02 from a 
set of 225 early-type galaxies in nearby clusters observed 
in the r-ba nd. This is consistent with the original obser- 
vation s by iDiorgovski fc Davisl lll987t) and IDressler et alJ 
Jl987h but noticeably different from the SDSS determi- 
nation (based on ~ 9000 early-type galaxies) of a — 
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1.49 ± 0.05 and b = -0.75 ± 0.01 feernardi et alJ l2003bl) . 
values more similar to the K-band fund amental plane of 
IPahre. de Carvalho. fc DiorgovskH l)l998f) . The reason for 
this discrepancy is not currently clear. 

Projections of the fundamental plane are also of interest 
in studies of galaxy evolution. The SDSS team measured a 
radius-luminosity relation of 



0.63±0.025 



Re OC L' 

and a Faber-Jackson (1976) relation of 



a oc L 



0.25±0.012 



(2) 



(3) 



Bernardi et alJ l2003al) . while work by Kauff mann et all 
J20031) on the ste llar mass content of SDSS galaxies allowed 
IShen et alJ l|2003f) to determine the scaling of effective radius 
with stellar mass: 



Re oc M" 



(4) 



The fundamental plane has typically been interpreted 
as a manifestation of the virial theorem, which relates a and 
R to the total enclosed (dynamical) mass Mdyn and predicts 
a correlation in the fundamental properties of galaxies: 

Mdyn I Mdvn \ ( L \ „ / Md 



2 

a oc 



R 



IR, 



giving the relation 

-i ^ Mdyn y 1 



Roc a I 



(•>) 



(6) 



Clearly this virial theorem expectation is incompatible with 
the observed fundamental plane in equation Q if light traces 
the dynamical mass (Mdyn oc L). A systematic variation of 
Mdyn i 'L with luminosity is usually assumed to "tilt" the 
virial theorem scalings into the observed ones. Since 

Mdyn I Md„n \ ( M* 



(7) 



this tilt could result from increasing M*/L with L due 
to, for instance, varying metallicity or stellar popula- 
tion age, or from increasing Md yn /M* with L due to 
higher dark matter fraction in the central parts of more 
luminous galaxies. Conflicting conclusions have been re- 
ported in the literature as to whic h of these explana- 
tions is correct (for instance, co mpare iGerhard et "ail 1200 ll 
with iPadmanabhan et alJ I2004F) . Structural or dynamical 
non-homology c ould also contribute to the fundamental 
plane tilt (e.g., |C ^aeJato^_d£_C^Jy^lhiL-j£-XlMl^£rj [l335j 
[Graham fc ColiessirH)97i r ip£Jire 

Numerical simulations are a powerful way to assess 
the dynamics of galaxy mergers. There is a substantial 
body of work devoted to simulations of mergers of galax- 
ies in the context of for ming ellipticals (for review s, see 
iBarnes fc Hernauisdll992l and iBurkert fc Naabl Eoolj) . Sev- 
eral studies have in particular addressed the question of 
whether collisionless mergers of galaxies preserve the funda- 
mental plane relations. These studies typically used simula- 
tions of ~ 10 4 particles (per initial galaxy) and have reached 
varying and sometimes conflic t ing r esults. For instance, 
iGonzalez-Garcia fc van Albadal i2003fl studied mergers of 
one-component systems (stellar bulges without dark matter 
haloes) with R e oc M ' 5 and found that the merger rem- 
nants ended up very near the fundamental plane defined 



by the progenitors. iDantas et alJ (|2003) found that mergers 
of bulge-only galaxies remained on the fundamental plane, 
while mergers of galaxies containing both bulge and dark 
matter components produced a fundamental plane with a 
tilt la rger than that observed. iNipoti. Londrillo. fc Ciottll 
(2003) used repeated mergers of one- and two-component 
galaxies to follow the merging process over a large range in 
mass. They found that neither the R e — M* nor the M* — <r e 
relation was preserved in their dissipationless mergers but 
the full fundamental plane of their merger remnants agreed 
with observations. 

In this paper we carry out dissipationless mergers of 
early type galaxies, focusing on major mergers of two- 
component (dark matter plus stellar bulge) systems. We 
thus do not address the origin of either bulge-dominated 
galaxies or the fundamental plane but rather investigate 
whether dissipationless mergers maintain the tight funda- 
mental plane relations, a scenario that is likely important 
for massive ellipticals. We pay particular attention to us- 
ing cosmologically motivated models that include adiabatic 
contraction of dark matter haloes and merger orbits drawn 
from cosmological simulations. We also use a significantly 
larger number of simulation particles than many previous 
studies, which is necessary to avoid spurious two-body re- 
laxation (and an accompanying artificial expansion of the 
stellar bulge). In the next section, we describe our initial 
models and choice of simulation parameters in more detail. 
Section |3] discusses properties of both the stellar and dark 
matter components of the merger remnants in our simula- 
tions, including their radial density profiles, triaxiality, if 
and how the fundamental plane scalings and structural ho- 
mology are preserved, and the changes in the mass ratio 
Mdyn I M* in the inner parts of the remnants. Section in- 
cludes a discussion and summary of our results. 



2 SIMULATIONS 

2.1 Galaxy Models and Initial Conditions 

Our model galaxies consist of two components, a stellar 
bulge and a dark matter halo. W e model th e initia l spatial 
distribution of the bulge with a iHernauistl il990l) density 
profile, 



p*(r) 



M, 1 



2na 3 r/a yl + r/a 



(8) 



with total stellar mass M„, quarter-mass radius r\u = a, 
and half-mass radius n/2 = (1 + \/2)a. Observationally, a 
useful radius is the effective radius R e , defined as the ra- 
dius of the isophote that encloses half of the total (pro- 
jected) stellar luminosity. For equation ®, the effective 
radius is related to the quarter- and half-mass radii by 
R e — 1.8153a = 0.752rx/2. The Hernquist profile in pro- 
jection follows the lde Vaucouleursl l)l948h r 1 ^ 4 la w to within 
35% in the radial range 0.06 < R/R e < 14.5 faernauisd 
1990), a region containing over 90% of the model's total lu- 
minosity. 

For the dark matter component, we consider 
two possible radial dens i ty pr ofiles: the NFW profile 
iNavarro. Frenk. fc White! Il997l) and a more realistic 
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version of this profil e that uses the a diaba tic contrac- 
tion approximation felumenthal et all Il98rj) to model 
the response of the dark matter halo to baryons. The 
uncompressed NFW profile is given by 

PcSc 



Pdm{r) = 



(cr/r v )(l + cr/r v ) 



(9) 



where the concentration parameter c is defined as the ratio 
of the virial radius r v to the radius at which the logarithmic 
slope of the density profile equals —2 (also known as the 
scale radius r s ). The additional parameter 8 C depends on c 
and the choice of virial overdensity At, relative to the critical 
density p c : 

A„c 3 



" c 3[ln(l + c) - c/(l + c)] ' (10) 

We use A,, = 200 throughout this paper. Since equation © 
is logarithmically divergent as r — > oo, we use the expo- 
nentia l truncation scheme described by ISprineel fc White! 
(f999) for our haloes at r > r v . Truncating the haloes in 
this manner requires « f 0% more simulation particles than 
a sharp truncation at the virial radius but keeps the outer 
part of the halo in equilibrium. 

We calculate the adiabatically compressed dark matter 
density profile from the initial (NFW) dark matter profile 
and the final baryon (Hernquist) profile under the assump- 
tions of homologous contraction, circular orbits for the par- 
ticles, and conservation of angular momentum r M{< r): 



[M b (n) + M dm {rP)] = r f [M b (r f ) + M dm {n)] . 



(11) 



Although it is not obvious that the adiabatic contraction ap- 
proximation should be directly applicable to elliptical galax- 
ies - it was originally proposed to describe the response of a 
dark matter halo to baryonic dissipation during disk forma- 
tion - simulations have shown it to be surprisingly accurate 
for mo deling ellipticals. Recent hydrodynamical simulations 
iGnedin et all20o3) find that equation 1111 overpredicts the 
amount of contraction at small radii but this difference is 
much smaller than the difference between haloes with and 
without adiabatic contraction, so we do not expect it to sig- 
nificantly affect our results. 

Fig- shows the initial density profiles (top panel) and 
1-d velocity dispersions (bottom panel) of the stellar bulge 
(symbols) and dark matter haloes (plain curves) in our sim- 
ulations. It also illustrates the effects of compression in the 
dark matter haloes (solid vs dashed curves). The top panel 
shows that the mass density interior to i? e = 2.8 kpc is dom- 
inated by the stellar component when the dark matter halo 
follows the NFW profile without compression, but compres- 
sion increases the central dark matter density significantly 
and raises it to be larger than the stellar density at r > 1 
kpc. Compression also steepens the NFW p a r" 1 density 
cusp to nearly isothermal (r~ 2 ) at r>l kpc (although the 
compressed density still asymptotes to r _1 at very small 
r). The bottom panel shows that adiabatic compression in- 
creases the velocity dispersion of both the dark matter and 
stellar components. It is also interesting to note that regard- 
less of compression, the mere presence of a stellar bulge in 
a dark matter halo increases the velocity dispersion of the 
dark matter component by up to a factor of ~ 2 (for r > 1 
kpc) in comparison with that of an NFW halo without a 
bulge (dotted curve). This effect also greatly reduces the ra- 



dius at which the "temperature inversion" characteristic of 
the NFW profile occurs. 

We construct the JV-Body models by setting the bulge 
and dark matter halo to be in equilibrium in the total 
gravitational potential. We assume that each component 
is both spherical and isotropic. The particle positions for 
each component are initialized from the density profile of 
that component. The particle velocities are drawn based on 
each component's equilibrium distribution function ,fi(E), 
which we calculate numer ically using Eddington's formula 
jBinnev fc Tremainelll987t) : 



ME) 



1 



V8tt 2 



d 2 pi dip 



dip 2 y/E-tp 



(12) 



where pi is the density profile of component i while ip 
is the total gravitational potential. (This is to be con- 
trasted with the typical one-component model, where p and 
ip are related by Poisson's equation.) This method of se- 
lecting the velocities is self-consistent and makes no as- 
sumptions about the local form of the velocity distribu- 
tion, allowing us to create stab l e initi al conditions (see 
iKazantzidis. Maeorrian. fc Moore) i2004T) for a full discus- 
sion of sampling velocities when constructing initial condi- 
tions for N-body models). Note that while the positions of 
the particles in each component are the same as would be 
given in a one-component model, the velocities are not due 
to the presence of the other component. 



2.2 Model Parameters and Merger Orbits 

Table^summarizes the parameters used in the galaxy mod- 
els and simulations for our series of production and test runs. 
Before adiabatic compression, the initial haloes in our sim- 
ulations have Mdm = 1O 12 M0 and concentration c = 10, 
giving a scale radius of r a = 16.26 kpc and virial radius of 
r v = 162.6 kpc. Gravity is the only physics in our simulation, 
meaning we can in principle interpret the results at different 
mass scales by a corresponding scaling of time and length. 
In practice, however, the presence of a stellar bulge com- 
plicates this scaling (see Section [I] for details). We consider 
two values for the total dark matter to stellar mass ratios, 
M dm /M, = 20 (M, = 5 x 10 10 M Q ) and M dm /M* = 10 
(M* = 10 ^Mq), which are cho sen to be representative o f 
the values estimat ed from SPSS jPadmanabhan et alJl20o"3) 
and weak lensing f§eliak 20o3). As discussed in Section |31 
we do not find the merger remnants to depend sensitively 
on this ratio. We relate the remaining free parameter of the 
Hernquist density profile - the length scale R e - to the stel- 
lar mass using th e observed R e — M * relation of early-type 
galaxies in SDSS <Shen et al.1120031) 1 : 



Re 



4.16 



M, 



10 11 M© 



kpc 



(13) 



Our choices of Md m /M* and R e do not guarantee that the 
velocity dispersion <r e of a model bulge will lie on the ob- 
served fundamental plane. In fact, observations of the fun- 
damental plane measure the Faber- Jackson L — o e relation 

Note the correct coefficient in Table 1 of lShen et alJ i2003t) is 
2.88 X 10~ 6 rather than 3.47 X 10~ 5 (S. Shen, private communi- 
cation) 
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Figure 1. Initial radial density profiles p(r) (top) and 1-d veloc- 
ity dispersion profiles o(r) (bottom) of the stellar bulge (symbols) 
and dark matter halo (plain curves) used in our simulations. The 
dark matter and bulge masses are set to be M^ m = 10 12 Mq 
and M» = here. An adiabatically compressed dark mat- 

ter halo (solid) has a larger p and a than the standard NFW 
halo (dashed) in the inner tens of kpc. The velocity profile of a 
standard bulge-less NFW halo of the same mass is plotted (dot- 
ted curve) for comparison. (The density profile of the bulge-less 
NFW halo is the same as the uncompressed NFW halo). 



while our simulations use the M* — cr e rela tion. We relate 
the tw o using the r* band values from SDSS feernardi et alJ 
l2003al) . Assuming M*/L = 3-3.5 in the r* band, consis- 
tent with the estimates of Kauffmann et al. (2003) , we find 
that a bulge of M* = 5 x 10 10 Mq should have a character- 
istic velocity dispersion of ~ 160 km/s. TableQshows that 
this is reasonably consistent with the values we obtain from 
model M20 (which includes adiabatic contraction) and is no- 
ticeably higher than the uncompressed halo model (M20u). 
In order to lie on the M* — cr e relation, the uncompressed 
model must have a very large value of M*/L ~ 8. Adia- 
batic contraction (and the associated increase in the dark 
matter fraction at small radii) therefore helps put the initial 
galaxy m odels closer to the observed fundam ental plane (see, 
however, iBorriello. Salucci. fc DaneseH2003l for an alternate 
view) . 

In addition to the internal structure of the bulge and 
dark matter halo, the orbit of a binary encounter is also 
important in determining the properties of the merger rem- 
nants. For instance, it is well-known that radial orbits tend 
to result in prolate remnants whereas orbits with significant 
angular momentum generall y give oblate or triaxial end- 
products iMoore et aljEooil and references therein). Most 
merger studies to date have assumed parabolic orbits with 
somewhat ad hoc impact parameters. We attempt to use 
more realistic orbits by drawing from recent results on the 
orbital distribution of merging dark matter haloes in cosmo- 
logi cal N-Body sim ulations. 

iBensonl (I2005T) . for example, has analyzed the orbital 
parameters of minor mergers in the Virgo Consortium's N- 
Body simulations of a ACDM model. 2 He finds that the 
orbital distribution peaks at se mi-major axis ~ r v and ec - 
centricity e slightly less than l. lKhochfar fc Burkertl (120051) 
studied similar quantities in major mergers (with mass ratio 
within 4:1) and found that the distribution of circularities 
peaks at e = yl — e 2 = 0.5, consisten t with the findings 
of Benson, a n earlier study bv | Tormenl dl997f> . and recent 
results from IZentner et alJ J2005I) . We therefore choose a 
"most probable" orbit defined by e = 0.5 and pericentric 
distance r v — 50 kpc, consistent with Khochfar & Burkert's 
relation r p = 210e 2 ' 07 kpc. We assign the radial and tangen- 
tial velocities v r and vg of this orbit to agree with Benson's 
analysis (his Figs. 2 & 7, and equation 2): v r = v v i r and 
vg — 0.7v vir where v 2 ir = GM v i r / Rvir- In addition to the 
"most probable" orbit, we also simulate a head-on parabolic 
merger for comparison. Finally, we perform one modified 
"most probable" simulation in which both v r and ve are re- 
duced by 15% in order to test how robust our results are to 
the choice of initial orbital velocities. The initial center-of- 
mass velocity for each run is listed in Tabled In all cases we 
start with the galaxy centers separated by two virial radii. 

2.3 Simulations and Resolution Studies 

We use GADGET dSpringel. Yoshida. fc Whitell20oih . an effi- 
cient parallelized tree code, to perform the simulations pre- 
sented below. We set the force softening e = 0.3 kpc for both 
star and dark matter particles; this choice corresponds to ei- 
ther 0.072i? e (M. = W U M @ ) or 0.11i? e (M* = 5x 10 10 M©) 

2 http:/ /www. mpa-garching.mpg.de/Virgo/ 
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Table 1. Simulation and model parameters of each galaxy in the 
initial conditions used in our production runs (upper 5) and test 
runs. In the run names, the merger orbits are distinguished by 
"M" and "R" for "most probable" vs. radial (parabolic) orbits; 
the mass fraction M^ m /M t is labeled by 20 vs. 10; the runs with 
initially uncompressed dark matter halos are denoted by "u" . 
Run N20 is identical to run M20 except the magnitude of each 
component of the initial orbital velocity is reduced by 15%. Nj, m 
and N* are the numbers of dark matter and stellar particles used 
in each simulation. All dark matter halos have virial masses of 
10 12 Mq and virial radii of 162.6 kpc. 



Run 



M . 



N, 



Ra 
e 



M10 


10 


5.5 X 10 5 


49000 


4.13 


173.5 


200 


M20 


20 


5.5 x 10 5 


24500 


2.85 


151.4 


200 


N20 


20 


5.5 x 10 5 


24500 


2.85 


151.4 


170 


M20u 


20 


5.5 x 10 s 


24500 


2.83 


122.7 


200 


R20 


20 


5.5 x 10 5 


24500 


2.83 


150.4 


250 


MlOtl 


10 


1.1 x 10 s 


49000 


4.13 


172.8 


200 


M20t0 


20 


1.1 x 10 5 


4900 


2.79 


151.3 


200 


M20tl 


20 


1.1 x 10 5 


24500 


2.76 


150.9 


200 


M20t2 


20 


1.1 x 10 5 


49000 


2.80 


150.1 


200 


M20tlu 


20 


1.1 x 10 5 


24500 


2.76 


123.1 


200 


R20tl 


20 


1.1 x 10 5 


24500 


2.76 


150.9 


250 



"effective radius of stellar bulge, in kpc 

6 aperture velocity dispersion of stellar bulge feci. 1151 . in km s _1 
c center-of-mass velocity for merger, in km s — 1 



10.0 




0.8 1.0 

(R/R.) 1/4 

Figure 2. Surface brightness profiles for run M20 at t=0 Gyr 
(dashed curve) and ten random projections at t=6 Gyr (solid 
curves) as a function of (R/ Re) 1 ^ and normalized to unity at 
R = R e (the initial profile is offset by 15% for clarity). A de 
Vaucouleurs law is a straight line on this plot. 



and stellar matter to resolve quantities of interest on scales 
larger than ~ 1 kpc. 



for our initial galaxy models. We run each model until the 
remnant reached virial equilibrium (between 4 and 6 Gyr de- 
pending on the orbit). Energy is conserved to within 0.5% 
for all runs and the global virial ratio 2T/\W\ is unity to 
within 1.5% at the final timestep in each case. 

We have performed several convergence tests to ensure 
that sufficient particle numbers and force resolution are used 
for both the dark matter and stellar components. We find 
that increasing the number of star particles A* from 24500 
to 49000 (runs M20tl and M20t2) at fixed dark matter par- 
ticle number does not affect our results for R e and <r e . It 
is, however, important to use Ndm ~ 5 x 10 s dark matter 
particles per halo since our test runs with Ad m ~ 10 5 show 
noticeable differences: for example, if we assume R e oc M", 
using Ad m = 10 5 rather than 5 x 10 5 causes an overesti- 
mation of a at the 15% level (for a compressed halo with 
Mdm/M* = 20). This effect can be understood in terms 
of two-body relaxation: using a small number of particles 
means that on small scales, individual particle interactions 
become important even when compared to the mean gravi- 
tational field. This causes an artificial heating, diluting the 
dark matter potential at the center of the model galaxy, and 
resulting in an artificially expanded stellar bulge. Similarly 
the final stellar velocity dispersions in the test runs tend to 
come out lower than in the production runs, again because of 
the shallower central potential. This effect leads to an over- 
estimate of R e and an under-estimate of a e , biasing the fun- 
damental plane correlations we hope to measure. As listed 
in Table Q our production runs all use 5.5 x 10° dark mat- 
ter particles and 24500 stellar particles (for Mdm/M* = 20). 
These values provide enough mass resolution in both dark 



3 PROPERTIES OF MERGER REMNANTS 

The global merger dynamics in our simulations using the 
"most probable" orbit is similar in all three cases (runs 
M10, M20, and M20u). As the galaxies spiral in toward each 
other, the dense stellar bulges stay intact while the more 
diffuse dark matter haloes mix rapidly. Dynamical friction 
then brings the bulges and the central core of the dark mat- 
ter haloes together and the merger is finished after approxi- 
mately 5 Gyr. The process for the radial orbit is about 30% 
faster even though the initial orbit is more unbound because 
the orbit brings the bulges together earlier. In the follow- 
ing sections we quantify the structure and kinematics of the 
merger remnants and then discuss the extent to which the 
merger preserves homology and maintains the fundamental 
plane relations of elliptical galaxies. 

3.1 Surface Brightness and Radial Density 
Profiles 

One distinctive property of the bulges in early-type galaxies 
is that their surface brightness profiles follow the de Vau- 
couleurs law, log I(R) oc -i? 1/4 . Fig. H shows I(R) for 10 
random projections of the final bulge from run M20 (solid 
lines) as a function of (R/Re) 1 / 4 . For comparison we also 
show I(R) of the initial bulge (dashed line; offset by 15% for 
clarity). Each of the random projections of the final bulge, 
as well as the initial profile, all follow the de Vaucouleurs 
law over the plotted range of e < R < Rqo (where Rqo is the 
projected radius containing 90% of the bulge mass). 
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Table 2. Fits to the radial density profiles of dark matter haloes and stellar bulges at t = and 6 Gyr in our five production runs. 
Parameters A, B, and C are defined in equation 1141 , Haloes are fit on the range [1,100] kpc, bulges on the range [0.32,16] kpc. 



Run 


time 


M 


Dark Matter Haloes 

A B 
(10 5 M Q kpc- 3 ) (kpc) 


L - 


M 


CH-ollar* Riil 

r l/4 

) ( k P C ) 






M10 





100.6 


14.8 


32.5 


1.8 


10.1 


2.24 


2.19 


0.97 




6 


133.8 


17.3 


34.1 


1.8 


18.6 


2.94 


4.47 


1.5 


M20 





100.4 


9.52 


39.4 


1.8 


5.02 


1.57 


1.55 


0.98 




6 


136.8 


12.1 


40.2 


1.8 


9.75 


2.26 


3.83 


1.6 


N20 





100.4 


9.52 


39.4 


1.8 


5.02 


1.57 


1.55 


0.98 




5.5 


144.1 


9.87 


44.8 


1.8 


9.82 


2.21 


4.31 


1.7 


M20u 





100.0 


134 


15.8 


0.98 


5.08 


1.56 


1.68 


1.1 




6 


136.2 


114 


18.8 


1.1 


9.92 


2.04 


3.88 


1.7 


R20 





100.4 


6.84 


45.2 


1.8 


5.06 


1.56 


1.55 


0.99 




4 


132.1 


14.6 


37.6 


1.7 


10.1 


2.89 


7.31 


1.9 



Although the remnants' projected density profiles 
mimic those of the initial conditions, the unprojected den- 
sity profiles are not necessarily identical. In order to quantify 
the changes in density structure induced by the merger, we 
fit both the bulges and dark matter haloes to the following 
profile: 

P(r) = {r/B) c (l + r/B)(- D - c >' (14) 

Since we find the outer slopes of both the dark matter haloes 
and stellar bulges to be mostly unaffected by the merger, 
we fix D to their initial values: D = 3 for the halo and 
D = 4 for the stellar bulge. M odels with D = 4 ha v e bee n 
investigated by iDehnenl jl99ot) and iTremaine et alJ lll994l) ; 
for these it is possible to relate A to B, C, and M t : A — 
(3 - C)M*/(47rB 3 ). Initially the bulges are all described by 
Hernquist density profiles, so A — M*/(27ra 3 ), B — a — 
r i/4i C — 1, and D = 4. The adiabatically compressed dark 
matter haloes differ from the original NFW profiles, so we 
also fit the initial compressed halo profiles to equation l|141 
with D = 3. 

The results are summarized in Table and show an 
interesting trend: the density profiles of the dark matter 
haloes retain their shape to a remarkable degree throughout 
the merger. In fact, the main difference between the initial 
and final pdm is an overall amplitude change by a factor of 
~ 1.4, regardless of the merger orbit or whether the halo is 
compressed. The bulges show more evolution after the merg- 
ers, especially in the inner slope of the density profile. In all 
cases the fitted inner slope of the bulge steepens noticeably. 
Although observations indicate t hat massive elliptical galax- 
ies have surface brightness cores IIFaber et al Jl997t) . the size 
of these cores tends to be ~ 100 pc (for the galaxy masses in 
our simulations). This scale is below our force softening of 
e = 0.3 kpc, so there is no inconsistency between our surface 
brightness or density profiles and observations of elliptical 
galaxies. As with the dark matter, the total profile for stars 
plus dark matter changes mostly by an overall amplitude 
and remains close to isothermal on the scale of R e . 

We also calculate the final virial radius of each dark 
matter halo and compare the enclosed mass M v with the 
initial value in Tabled In each case only about 65-70% of 
the initial virial mass ends up in the final virial radius; a sim- 
ilar ratio was found in our earlier study of dark-matter-only 



mergers feovlan-Kolchin fc Mall2004l) . This "puffing up" of 
the outer part of the dark matter halo is in stark contrast 
to the stellar bulge, where essentially all of the initial stellar 
mass M* ends up in the final bulge. (This result is insen- 
sitive to changes in merger orbit, compression of halo, and 
the value of Mdmj 'M».) The ratio of the dark matter to stel- 
lar mass within the virial radius therefore decreases after a 
major collisionless merger. However, as discussed below in 
Section r3.5l the dark matter-to-stellar mass ratio within the 
central region of a galaxy (within R e ) shows the opposite 
behavior. 

3.2 Remnant Shapes 

We find that the merger remnants in our simulations tend 
to be triaxial, with their structure determined by the initial 
orbit: radial orbits produce nearly prolate remnants while 
orbits with non-zero angular momentum yield triaxial prod- 
ucts. 

To quantify this trend, we determine the axis ratios of 
the final products using an iterative eigenvalue technique 
jDubinski fc Carlberdll99lT) . We diagonalize a modified in- 
ertia tensor Iij = ^ XiXj/R 2 (where R 2 = x 2 + y 2 /(b/a) 2 + 
z 2 /(c/a) 2 is the ellipsoidal radius of a particle and a> b > c 
are the eigenvalues of Iij) by first assuming spherical symme- 
try and computing Iij , then using the eigenvalues of iteration 
n as input values for iteration n + 1. We also compare the 
axis ratios calculated both cumulatively {r\ < R) and in dif- 
ferential shells (Ri < Ti < 7?2). For the dark matter haloes, 
these two method give similar results at most radii, but near 
r v the cumulative method is less sensitive to local changes in 
the axis ratios due to the contribution of many particles at 
smalle r radii, an effect noted previously in lKazantzidis et alJ 
(2004). The cumulative method is, however, more robust for 
the bulges due to the relatively low number of stellar parti- 
cles. For the following analysis we thus use the cumulative 
method for measuring the shapes of both components. 

The top two panels of Figure [!|] show the axis ratios of 
the remnant stellar bulges at t = 6 Gyr from four of our 
production runs (we do not include run N20 because its axis 
ratios are nearly indistinguishable from run M20). The ra- 
dial orbit (dotted curves) results in a prolate bulge with 
both b/a and c/a m 0.5, irrespective of radius. The "most 
probable" orbit, on the other hand, yields a triaxial bulge 
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Figure 3. Triaxiality of the merged stellar bulges (top) and dark 
matter haloes (bottom) as a function of radius at t = 6 Gyr from 
four production runs. Axis ratios b/a (left) and c/a (right) are 
shown (a > b > c by definition). Radial orbits (dotted curves) 
tend to produce prolate bulges and haloes, whereas orbits with 
non-zero angular momentum result in triaxial bulges and nearly 
oblate haloes. Note that the ratios for the bulge are plotted lin- 
early in radius, while those for the halo are plotted logarithmi- 
cally. 



with axis ratios that have a mild radial dependence. These 
axis ratios are in agreement with a recent analysis of the 
SDSS data that found that bulges of early-type galaxies ar e 
consistent with tria xial ellipsoid s dVinc ent fc RvderJl2005l) . 
though we note that lAlam fc Rvde rj(20fj) also found that 
prolate ellipsoids are consistent with the SDSS early data re- 
lease sample. Simulations of a brightest c luster galaxy from 
a cosmological simulation |Dubhrskill99ct) find a stellar rem- 
nant with b/a = 0.66, which is similar to our remnants of the 
"most probable" orbits, and c/a = 0.47, which is closer to 
what we find for a radial merger. The latter is a consequence 
of the correlated directions of Dubinski's multiple mergers 
along filaments. 

The bottom two panels of Figure [3] show the axis ra- 
tios for the merged dark matter haloes. The radial orbit 
again results in a mostly prolate halo, whereas the three 
"most probable" orbit runs give nearly oblate haloes with 
b ~ a and c/a ~ 0.7, largely independent of Mdm/M, 
and halo compression. It is useful to compare our results 
w ith those from dark-matte r-only cosmological simulations 
of lBailin fc Steinmetzl i2004) . For a fair comparison, we cal- 
culate the axis ratios of our remnants in a shell from 0.25 
to 0.4 r v , the same range as Bailin and Steinmetz use. We 
find the "most probable" orbit to give b/a ~ 0.9 — 0.95 
and c/a ~ 0.7 while the radial orbit gives a r emnant with 
b/a ~ c/a ~ 0.7. Both agree with the results of ljing fc Sutol 
(2002) and Bailin & Steinmetz (cf their fig. 5), though our 
values of b/a tend to be at the high end of their distribution. 

The difference in the "most probable" orbit runs be- 
tween the bulge (triaxial) and halo (oblate) shapes is note- 



worthy, as it relates to the difference between the merging 
of the haloes and of the bulges. The haloes merge relatively 
quickly, but the bulges circle each other and spiral inward 
due to dynamical friction on the haloes. The final plunge of 
the bulge tends to be relatively radial even though the initial 
orbit has substantial angular momentum, giving the rem- 
nant bulge a different shape from the remnant haloes. Quan- 
titatively, we find that the amount of (initial orbital) angular 
momentum lost from the bulge is ~ 75% (i.e. Lf m 0.25Li). 
As a related measure, we calculated the ratio of rotation ve- 
locity to velocity dispersion for several projections. As might 
be expected, the maximal v rot /a ~ 0.5 at R s occurs in the 
plane of the orbit, while for most lines of sight v rot /o < 0.3, 
typical of a slowly rotating, pressure-supp orted bulge. This 
is con sistent with recent simulations by iNaab fc Burker tl 
(2003), who found that equal mass mergers of disk galaxies 
tend to result in pressure-supported ellipticals while merg- 
ers with mass ratios of 3:1 and 4:1 result in oblate rotators 
(with v ro t/a > 0.7). 



3.3 The Fundamental Plane Relations 

Before investigating the scaling relations among the effec- 
tive radius R e , velocity dispersion a e , and mass of the stel- 
lar bulge M*, we first describe how we determine R e and 
cr e - Because the remnant bulges are triaxial, their observ- 
able properties depend on the angle from which the bulge is 
viewed. To quantify this effect, we calculate the properties of 
each of our remnant bulges for 10 4 random viewing angles. 
A standard technique for obtaining the effective radius is to 
fit the projected profile to a Sersic profile and use the de- 
rived R e as the half-light radius. As described in Appendix 
lAl however, we find this procedure does not give unique re- 
sults, so we employ an alternate technique here. We calculate 
the surface mass density perpendicular to each line of sight 
in spherical bins as well as the projected enclosed mass at 
each radius. We define the projected radius of the particle 
enclosing half the total mass to be the half-light radius R e . 
This procedure recovers the correct R s for the initial Hern- 
quist profile to within 1-2 %. For the velocity dispersion, we 
use the surface- brightness weighted dispersion within R e , 



rR, 

e\ = o 2 a {R e ) = ifi- 



„2 



(R)I(R)RdR 



f* e I(R)RdR 



(15) 



where I(R) is the surface brightness and <r; os is the line-of- 
sight velocity dispersion. Observationally, a commonly used 
definition is the dispersion measured within an aperture of 
size R e / 8. For our simulations, however, R e /8 is compara- 
ble to the force softening e, so we instead define our veloc- 
ity dispersion within R e . We do not expect this modifica- 
tion to bias our results with respect to observations because 
the apertu re dispersion profile a a (R) is generally quite flat 
within R e jGerhard et alJl200ll) . Any error introduced will 
almost certainly be less than that associated with empiri- 
cally correcting the observed aperture dispersions to R e /8. 
To be conservative with respect to force accuracy, we use 2e 
as the minimum radius for the integral in equation 1151 . 

Figure [I] shows the resulting distributions of Re (top 
panels) and a e (middle panels) from 10 4 random view- 
ing angles of our merger remnants from two simulations 
with different orbits: "most probable" (left panels) and radial 
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Table 3. Mean (left) and mode (right) of the distributions of the effective radius R e and velocity dispersion <r e from 10000 random 
projections of the remnant bulge in each of our five production runs. Parameters a and j3 are the resulting exponents for two projections 
of the fundamental plane relation: R e oc M" and M* oc erf. 







Mean 








Mode 






Run 


Re 




a 


& 


Re 




a 






(kpc) 


(km/s) 






(kpc) 


(km/s) 






M10 


6.94 


200.0 


0.75 


4.9 


6.45 


204.5 


0.64 


4.2 


M20 


4.80 


175.7 


0.75 


4.7 


4.58 


179.0 


0.69 


1.1 


N20 


4.85 


173.0 


0.77 


5.2 


4.63 


176.0 


0.70 


1.6 


M20u 


4.53 


153.0 


0.68 


3.1 


4.49 


150.0 


0.67 


3.5 


R20 


6.63 


162.8 


1.23 


8.7 


7.01 


154.3 


1.31 


27 
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Figure 4. Distributions of R e (top) and <r e (middle), as well as 
the correlation at fixed mass between Re and cr e (bottom), for 
10000 random viewings of the merged stellar bulge in runs with 
Mdm — 20M*. Left: "most probable" orbit with vg/v r ~ 0.7 (run 
M20). Right: radial orbit (run R20). Note that the mean (dotted 
line) and mode can differ significantly, as the distributions are 
heavily skewed. The distributions are also dependent on orbit: 
for run M20 the mean R e is greater than the mode and the mean 
a e is smaller than the mode, while the reverse holds for run R20. 
In both cases the "observed" R e and er e are highly correlated, as 
indicated by the bottom two panels. 



(right panels). The distributions are highly non-Gaussian, 
although in each case there is a reasonably well-defined mode 
that differs significantly from the mean (indicated by the 
vertical dotted line). The two orbits result in strikingly dif- 
ferent distributions of R e and <r e . This difference is primarily 
because the two classes of orbits yield remnants with differ- 



ent shapes (see Section 13.21 above). The "most probable" 
orbit yields a triaxial remnant for which one is more likely 
to observe a small R e whereas the radial orbit yields a pro- 
late remnant for which one is more likely to observe a large 
Re. In addition, mergers on orbits with non-zero angular 
momentum (such as the "most probable" orbit) suffer more 
dynamical friction and thus more energy transfer from the 
stellar bulge to the dark matter (see also Section s-inl and Ta- 
ble HJ. This effect leads to a smaller and more tightly bound 
bulge with a larger velocity dispersion. 

The bottom panels of Figure^Jshow the correlation be- 
tween the fundamental plane parameters Re and cr e for fixed 
stellar mass from the 10 4 viewing angles. It is evident that 
the measurements are highly correlated; viewing angles that 
yield a relatively small effective radius result in a relatively 
large velocity dispersion. This effect is due to projection: 
looking down the major axis of a prolate remnant, for ex- 
ample, leads to a smaller effective radius and larger velocity 
dispersion than if the viewing angle is along a minor axis. 
These correlations are also expected from the fundamental 
plane at fixed luminosity (or mass): from the SDSS values 
for equation R e oc <r~ 3 at fixed luminosity. Although 
the fundamental plane cuts through the R e — er e plane in 
a non-trivial manner, making a direct comparison between 
the spread in our R e — a e correlations and the scatter of the 
fundamental plane difficult, the thickness of the correlations 
we present here are consistent with the observed small scat- 
ter in the fundamental plane. However, it is not clear how 
additional effects, such as scatter in the global Mdm/M t and 
in initial sizes and shapes of merging galaxies, would affect 
this scatter. 

Using the results of Figure [I] we can now assess the 
extent to which the stellar bulges remain on the projections 
of the fundamental plane after major dissipationless merg- 
ers. Table El lists the mean and mode of R e and cr e for the 
merged stellar bulges and the slopes for power-law fits to 
R e — M* and M* — cr e from our five production runs. Since 
the mean and modes differ significantly, we present results 
for each definition of the "average" remnant. The results for 
R e - M* shown in Table 01 are more robust than those for 
M* — (J e because of the steepness of the latter relation. 

Table |31 shows that for the "most probable" orbit 
with adiabatically compressed dark matter haloes (i.e. runs 
"M10" and "M20", which are the best motivated models 
studied in this paper), the derived fundamental plane slopes 
using the mode are reasonably consistent with observations: 

Re oc M" , a = 0.64 - 0.69 , 
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Table 4. Parameters from the bulge energy conservation equa- 
tion fean. USt derived from the simulations. We quote numbers for 
the remnant bulge structural parameter ft using both the mean 
R e and mode R e obtained from 10 4 random projections (listed 
in Table [3) . Parameters f oro and ft characterize the initial or- 
bital energy (at virial radius) and the energy transfer between 
the stellar and dark matter components due to mergers. 



Run 


fi 


ff,rnean 


ff,mode 


forb 


/* 


M10 


0.342 


0.385 


0.358 


0.118 


0.344 


M20 


0.367 


0.444 


0.424 


0.144 


0.469 


N20 


0.367 


0.448 


0.428 


0.215 


0.394 


M20u 


0.265 


0.319 


0.316 


0.144 


0.375 


R20 


0.367 


0.489 


0.517 


0.00715 


0.183 






OC <7 e , 


f3 = 4.1 


-4.2. 







(16) 

Recall that the initial values of R e , M*, and a e for these 
simulations also lie on the fundamental plane (see Section 
12.21 . In comparison, the "most probable orbit" model with 
an uncompressed dark matter halo (run "M20u") produces 
a M* — <7 e relation somewhat shallower than the observed 
relation, while radial orbits (run "R20" ) yield a remnant well 
off the fundamental plane (/3 ~ 27!). In fact, the remnant 
from the head-on merger is much further off the fundamental 
plane than the intrinsic scatter of the plane, suggesting that 
few major mergers are radial or nearly so. This inference is 
consistent with the low frequency of head-on mergers in the 
orbital analyses discussed in Section f2.2l 



3.4 Homology and Energy Conservation 

Semi-analytic galaxy formation models often use energy con- 
servation arguments and rely on the assumption of homology 
for predicting how elli ptical galaxies evolve durin g dissipa- 
tionless mergers <Cole et all200d:IShen et al.l2003l) . The en- 
ergy conservation equation is usually derived using the fact 
that on dimensional grounds we can quantify the energy of 
a single bulge of mass M* as 



E = -f 



GMj 
R* 



(17) 



for some characteristic radius R, and structural parameter 
/ that depends on the properties of the galaxy model under 
consideration. For a bulge with the Hernquist profile (with- 
out a dark matter halo), E = -GM*/12a, so / = R e /12a = 
0.151 using 7?, = R e , or / = (1 + \/2)/12 = 0.201 using 
R t — Tx/2- For our various bulge- plus-halo models, we find 
/ = 0.265 — 0.367 (see Table[4j due to the presence of a dark 
matter halo. 

Using equation 1171 . we can write an energy conser- 
vation equation for the bulges in a binary merger of two 
galaxies (with initial bulge masses Mi and M2) as 



Rf 



M{ ..Mi MrM 2 

h lh +h lh + {Ub + ft) -R7Tlh ' (18) 



where f oro contains information about the orbital energy 
and ft describes the energy transfer between the stellar and 
dark matter components. With this definition, ft > implies 
that the final bulge is more tightly bound than the initial 
bulges, i.e. energy is transferred from the bulge to the dark 



matter halo. (An equation analogous to equation 1181 can 
be used to describe the dark matter component.) As noted 
by previous authors, setting ff — f\ = fi and f or b = ft = 
(i.e., perfect homology and a parabolic orbit with no energy 
transfer between the dark matter and the stellar component) 
implies R oc M* for an equal mass merger, which does not 
match the observed relation. Alternatively, assuming homol- 
ogy is satisfied (// = fi) and R oc M", one can derive that 



forb + ft 

fi 



r + 1 



[(C + i) 



-r 



-1 



(19) 



for binary mergers between galaxies with mass ratio £. The 
observed R e — M* relation is a = 0.56, implying f orb + ft = 
1.43/i for equal-mass mergers if homology is satisfied. 

To assess homology in our simulations we calculate the 
values of / using the energy defined in Appendix ITU 



E 



E 



E 



Grriiinj 



— Gm,i 



(20) 



where m = \ fi — rj-| and i and j run over all stellar particles 
while k runs over all dark matter particles. The resulting 
values for / are shown in Table [I] Homology is preserved at 
the ~ 15% level in all three runs with the "most probable 
orbit" (runs M10, M20, and M20u): / increases slightly in 
these runs, with the mean value changing more than the 
mode, but in each run both values are within 20% of the 
initial value. The same is true for run N20, which is identical 
to M20 aside from having lower orbital velocity. The radial 
merger (run R20), however, breaks homology significantly 
more than the "most probable" orbit. In all five production 
runs, the ratio of bulge-dark matter to bulge-bulge potential 
energy increases, reflecting the relative distribution of dark 
matter and stars, a point which we discuss further in Section 
1331 

Table 0] also gives the values of f or t and ft, the pa- 
rameters characterizing the amount of initial orbital energy 
(at the virial radius) and the amount of energy transferred 
between the stellar and dark components. For the "most 
probable" orbit, f t ~ 0.34 — 0.47, indicating that there is 
a significant amount of energy lost from the stellar bulge 
to the dark matter halo during the merger (partially via dy- 
namical friction). The energy transfer is significantly less for 
the head-on merger in which dynamical friction is less im- 
portant. The values of (f t + forb)/ fi ~ 1.3 — 1.7 for the most 
probable orbit are reasonably consistent with the estimate 
above for the amount of energy transfer required to explain 
the observed R e — M* relation. Interestingly, a comparison 
of run N20 to M20 show that while N20 has less orbital en- 
ergy, it also less energy transferred from the bulge to halo 
and the net result is a halo with nearly identical homology 
properties to those from run M20. 

3.5 The Dark Matter Fraction 

As noted in Section the fundamental plane relations are 
tilted from the simple one-component predictions based on 
the virial theorem. One possible explanation for the origin 
of this tilt is a luminosity-dependent dark matter fraction 
within the effective radius. In this section we test this hy- 
pothesis by quantifying how the dark matter fraction evolves 
during our simulations. 
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Figure 5. Ratio of dark matter mass to total mass interior to 
radius r as a function of r for initial condition (solid curve) and 
remnant system (dashed curve) for simulation M20 (compressed 
halo; upper curves), M10 (compressed halo; middle curves), and 
M20u (uncompressed halo; lower curves). The symbols mark the 
effective radius R e for the given simulation. In each case the frac- 
tion of dark matter interior to R e increases after the merger. 



In Figure |K| we plot the ratio of dark matter mass to 
total mass Md m (< r) / M to t(< r) as a function of radius for 
the three simulations using the "most probable" orbit: runs 
M20 (upper curves) and M10 (middle curves) have com- 
pressed dark matter halos while run M20u (lower curves) 
has a standard NFW halo. Both the initial (solid curves) 
and remnant (dashed curves) galaxies are shown and the 
mode effective radius R e of each model's bulge is marked 
with an asterisk. The results for the initial models (solid 
lines) show the dramatic effect of adiabatic contraction in 
the inner region of a galaxy: the dark matter mass inte- 
rior to R e in the adiabatically contracted model is approxi- 
mately twice that of an uncompressed NFW halo. After the 
merger the dark matter fraction interior to a given radius 
decreases, but the dark matter fraction interior to the final 
Re is greater than the dark matter fraction interior to the 
initial R e for both simulations. This increasing dark matter 
fraction with stellar mass is consistent with that needed to 
explain the observed tilt in the fundamental plane. Quanti- 
tatively, we find Md yn oc Ml' 11 for the compressed haloes 
and Mdyn oc Ml 19 for the uncompressed haloes. These re- 
su lts are reasonably con s istent with the observational results 
of iPadmanabhan et alJ i2004) from SDSS. The amount of 
dark matter within R e can be constrained in other ways, 
particularly through gravitational lensing. Though a large, 
statistically significant sample of lenses is needed, the dark 
matter fraction we obtain is co nsistent with analyses of 
existing lensing observat ions (e.g. iTreu fc Koopmanal2004t 
iRusin fc Kochanekll2005l) . 

One way to understand the increasing dark matter frac- 
tion is based on the dark matter mass profile in the in- 
ner parts of the galaxy: assume that the dark matter mass 



profile is given by Mdm(r) oc M v r v . For an NFW profile, 
r\ « 2 at the small radii appropriate near R e . For an adi- 
abatically compressed halo, r\ w 1.2 (Tabled- If the bulge 
stays roughly on the R e — M* relation, R e increases by a fac- 
tor of « 1.5 during an equal-mass merger. The stellar mass 
within R e increases by a factor of 2 during such a merger, 
while the dark matter mass within R e increases by a factor 
of ~ 2 * 1.5' 7 ~ 3.25 — 4.5 if the virial mass also doubles. If 
the final virial mass is only a factor of « 1.3 of the initial 
virial mass (Table |5J, the dark matter mass within R e in- 
creases by a smaller factor of w 1.3* 1.5 V w 2.1 — 3. All of 
these estimates yield a dark matter fraction within R e that 
increases after the merger, consistent with our simulations. 



4 SUMMARY & DISCUSSION 

We have performed numerical simulations of dissipationless 
mergers of equal-mass galaxies containing a stellar bulge and 
a dark matter halo to study how observable properties of 
elliptical galaxies scale during mergers. Some key features 
of our simulations include: 

• We consider both standard NFW dark matter haloes 
and NFW haloes modified by adiabatic response to baryonic 
dissipation. 

• We set up the initial stellar bulges using the observed 
Re — M* projection of the fundamental plane (most of which 
also follow the observed M* — a e relation). 

• Our primary merger orbit is taken from the most prob- 
able value of a cosmological distribution, which is weakly 
bound with ve m 0.7 v r . We also simulate a merger on a ra- 
dial parabolic orbit and a merger on an orbit with the initial 
orbital velocities reduced by 15% for comparison. 

• We perform a large number of random projections of 
each merger remnant in order to determine the effects of 
viewing angle on the observable characteristics. 

• Our resolution tests indicate that ~ 5 x 10 5 dark matter 
particles per halo are necessary to accurately resolve the dy- 
namics of the stellar component on the scale corresponding 
to observations. Test runs with a smaller number of particles 
lead to artificially "puffed up" stellar remnants. 

The main results of our work are: 

• As in previous studies, we find that the remnant bulges 
retain a de Vaucouleurs profile. In addition, the spherically- 
averaged radial density profiles of the dark matter haloes 
appear to be scaled-up versions of those of their progenitors, 
i.e., with the dark matter density increasing by a simple 
normalization factor. 

• The remnant shapes depend strongly on the orbit of the 
merger. The cosmologically "most probable" orbit results in 
a nearly oblate halo with a triaxial bulge, while less probable 
radial orbits result in prolate bulges and haloes. 

• The Re — M* and M» — er e projections of the fundamen- 
tal plane are roughly preserved for an "average" viewing an- 
gle in major dissipationless mergers with the most probable 
orbit, independent of the global ratio of dark to luminous 
matter (Table 0. There are, however, significant variations 
with viewing angle, and the observables are strongly corre- 
lated (Fig. The viewing-angle scatter appears consistent 
with the observed scatter in the fundamental plane. Our 
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simulations thus suggest that the fundamental plane is con- 
sistent with models in which massive elliptical galaxies un- 
dergo significant growth via late dissipationless mergers. If 
the merging elliptical galaxies contain central massive black 
holes, and if black hole masses roughly add when they co- 
alesce (i.e., there is not significant energy loss to gravita- 
tional waves), then our results also i mply that dissipation- 
less mergers maintain the Mbh — cr e iGebhardt et al. 2000; 
Ferrarese & Mcrritt 2000) and M BH - M* iMagoTrianet^dl 
ll99Sl:lHaring fc Rixll200*i relations. 

• The merger remnants contain more dark matter (rel- 
ative to stellar matter) within R e than their progenitors: 
Mdyn oc M*' 11-1 ' 19 , consistent with the hypothesis that the 
tilt of the fundamental plane is due to an increasing dark 
matter fraction within R e . 

• An appreciable amount of energy and angular momen- 
tum are transfered from the stars to the dark matter in the 
merging process. Quantitatively, the energy transfer found 
in our simulations is similar to that needed in semi-analytic 
galaxy formation models to explain the fundamental plane. 

One important goal of studies of elliptical galaxies is to 
understand the origin of the tilt in the fundamental plane. 
The virial theorem expectation relates the stellar veloc- 
ity dispersion and effective radius to the total mass within 
R e : (Je oc Md V n(< Re) / Re- The fundamental plane relates 
stellar velocity dispersion, luminosity, and effective radius: 
R e oc a~ 3 L 3 / 2 . As discussed below equation (J7J, if metal- 
licity or stellar age differences are solely responsible for the 
fundamental plane tilt then Md yn oc M* and M» oc L x with 
x w 1.15. On the other hand, if dark matter fraction is re- 
sponsible, then A/* oc L and Mdyn oc M* with x w 1.15. 
Our simulations contain no star formation or metallicity in- 
formation and thus M* oc L by assumption; we therefore do 
not have a direct way to evaluate the importance of metal- 
licity variations (thoug h we note that there is observational 
evide nce for them; e.g. iTrager et alj|200ot ICappellari et al] 
120051) . Nevertheless, we find that because of a change in the 
amount of dark matter with R e due to a merger, dissipation- 
less (major) mergers move elliptical galaxies roughly along 
the fundamental plane. It is interesting to note that the 
nearly homologous nature of the bulges in our simulations 
indicates that our results may hold for unequal mass merg- 
ers as well: using mass ratio £ = | rather than £ = 1 results 
in only a 1% change in (f or b + ft)/ fi from equation I I19I . 

Because our simulations involve only gravitational in- 
teractions, the results can be interpreted at different mass 
and length scales by an appropriate rescaling of time. This is 
somewhat complicated in two-component simulations such 
as ours, however, because bulges and dark matter haloes 
do not scale in the same manner: for example, R e oc M» ' 56 
while r v oc M° 33 . If we rescale the length and mass ac- 
cording to the virial theorem definition, then the bulge no 
longer lies on the fundamental plane. For example, con- 
sider a mass scale 10 times larger than what is presented 
in this work (keeping the concentration constant). Both the 
virial radius and effective radius increase by a factor of 
10 1/3 ; then for run M20, the bulge effective radius will be 
R e = 10 1/3 x 2.81 kpc = 6.05 kpc rather than the 10.24 
kpc predicted by the SDSS scaling (equation This dif- 
ference is comparable to the 1 — a scatter Shen et al. find 
in the R e — M* relation. As a result, we can indeed scale 



our simulations up or down by a factor of 10 in mass and 
10 1 ' 3 in radius and still be near the fundamental plane. Since 
the rescaling only affects the amplitude of the relations, the 
slopes in equation 1161 and Tableware actually preserved 
and the scaled mergers would move along the fundamental 
plane but with a constant offset. 

The results of Section E-i.2l show that our remnant bulges 
are either triaxial (non-radial orbit) or prolate (radial orbit), 
in agreement with analyses of SDSS data. In this context it is 
important to point out that our initial conditions are spher- 
ical and isotropic, while both the remnants of our mergers 
and observed ellipticals tend to be triaxial. It would cer- 
tainly be desirable to use more realistic (i.e. non-spherical) 
initial conditions in future work, though there are signifi- 
cant difficulties in doing so: operationally, it is not possible 
to create stable triaxial initial conditions corresponding to a 
desired density profile. It would also be of interest to extend 
the results presented here by including a gaseous or stel- 
lar disk. Such simultions should help shed light on the ori- 
gins of both spheriodal systems and the fundamental plane, 
topics that we have not addressed in this work. We note 
that in simul ations o f merg ers of halo-bulge-disk systems 
performed by iBarnesI lll992l) , the stellar remnant followed 
a similar radius-mass relation (r oc M ' 5 ) to what we find 
for halo-bulge mergers, perhaps implying that inclusion of a 
collisionless stellar disk would not change the conclusions of 
this work drastically, but this issue must be resolved with 
future simulations. 
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APPENDIX A: SERSIC FITTING 

A common paramete rization of the surface brightness profile 
of a galaxy is that of Sersic ( 1968): 

I(R) = I(R e ) e X p{-b(m)[(R/R e ) 1/m - 1]}. (Al) 

Though the function b(m) cannot be expressed in 
clos ed form, an ac c urate asymptotic expansion is given 
by ICiotti fc Bertinl (ll99St) as 6(m) w 2m - 1/3 + 
4/(405m); bulges and early- type galaxies ar e typically well- 
approximated by the Ide Vaucouleursl 1^948) law of m = 4. 
Fitting an observed surface brightness distribution to a Ser- 
sic profile is a standard technique for determining the (pro- 
jected) half-mass radius R e of a galaxy. We find that fitting 
our projected density profiles to Sersic profiles is not ro- 
bust, however: changing the radial range of the fit leads to 
systematic variations in both the Sersic index and the cor- 
responding effective radius. 

This dependence of fitting results on the range of radii 
can be seen for the initial bulge that is assigned a Hernquist 
profile. Fig. IA1I shows the derived R e (top panel) and m 
(bottom panel) for an input Hernquist profile with R e = 2.81 
kpc as a function of the maximum radius of the fit. Two 
different values for the minimum radius of the fit are plotted: 
0.1 kpc (X symbols) and 0.35 kpc (square symbols). In both 
cases the recovered R e and m are increasing functions of 
the maximum radius of the fit; additionally, the recovered 
R e is lower than the input R e = 2.81 kpc and the Sersic 
index m is less than 4 for all reasonable maximum values of 
the fit. This result is not unexpected since Hernquist models 
deviate from the de Vaucouleurs profile at both small and 
large radii, but it does present a problem for quantitatively 
characterizing the stellar bulges in our simulations since the 
fitted R e for the initial Hernquist profile can change by 20% 
or more depending on the range of radii considered; the final 
bulges show similar trends. In fact, Blanton et al. (2005) cite 
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Figure Al. A demonstration of the dependence of R e and m 
obtained from Sersic fitting (equation IA1I on the radial range 
of the fit for a Hernquist profile with R e = 2.81 kpc. The top 
panel plots the derived effective radius R e as a function of the 
maximum radius to which the fit is performed while the bottom 
panel plots the fitted Sersic index m. The best-fit values of both 
R e and m depend sensitively on both the maximum (horizontal 
axis) and minimum (symbols) radii used in the fit. 



a similar issue with Sersic fits to SDSS observations (see 
their appendix A and fig. 10). 

As a result, we do not use Sersic fits to determine the 
effective radii of either our initial conditions or merged el- 
lipsoids. Instead, we take a line of sight and calculate the 
surface mass density perpendicular to the line of sight in 
spherical bins. We also calculate the projected enclosed mass 
at each radius. We define the projected radius of the parti- 
cle enclosing half the total mass to be the half-light radius 



(R e )- This procedure recovers the correct R e for the initial 
hernquist profile to within 1-2 %. 



APPENDIX B: ENERGY CONSERVATION 

Consider a system of particles with distribution function / 
obeying the collisionless Boltzmann equation 



at 



+ v 



df 

dx 



dx 



df 

dv 



= 0. 



(Bl) 



By taking, in turn, a velocity and spatial moment of 
equation IB II we can ob tain the tensor virial theorem 
l|Binnev fc Tremainelll987i sec. 4.3): 



'a/3 



2 dt 2 
where 
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a/3 • 



a/3, 



K, 
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w a0 



pXaXpd X, 



~ / pv ci v l3 d 3 x, 



d$ 3 
pxfj- — d x 

OX a 
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(B3) 
(B4) 
(B5) 



and v a Vj3 is the mean of v a vp. For a steady-state system, 
the second derivative of the inertia tensor I a g vanishes and 
the trace of equation <B2t gives the standard scalar virial 
theorem 2K — —W. If we consider a system with a bulge 
(b) and dark matter (d), each with a distribution function 
obeying the collisionless Boltzmann equation with respect 
to the combined potential 3 , we can write the same tensor 
virial theorem for each component by replacing the total 
mass density p with the mass density of that component. 
This modification is straightforward for the inertia tensor 
and kinetic term but is somewhat more subtle for the po- 
tential term W a p- By writing the potential as the sum of 
the potential of each subsystem, 4> = + "Iv, the trace for 
the bulge breaks into two separate parts: 



W b 



- d® b 3 
PbX ■ -^-d x - 
ox 



PbX ■ 



dx 



d ,3 



d?x. 



(B6) 



The first term in the above equation, which we denote Wbb, 
is the familiar self-energy of the bulge. The second term, 
Wbd, describes the interaction energy of the bulge with the 
dark matter halo. If we write the mass density as a sum over 
(Dirac) delta functions 



Pb 



m,iS(x - 



(B7) 



and similarly for component d, then this term takes a more 
transparent form: 



(Si-Xj) 



(B8) 



This quantity Wbd is energy belonging to the bulge due 
to the interaction of the bulge with the dark matter such 

3 Since phase space densities are additive, the total phase space 
density ftot = fb + fd a l so obeys the collisionless Boltzmann 
equation. 
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that the total potential energy of the bulge is defined as 
Wt = Wbb + Wbd- [The second term in equation 1201 is Wbb 
while the third term is Wbd-] We similarly define the total 
potential energy of the dark matter halo to be Wd = Wdd + 
Wdb- Note with this definition we have 



which is simply the "cross" potential energy of the bulge 
and dark matter components, and thus 



where Wtot is the total potential energy of the system. For 
each galaxy used in our initial conditions, 2 Tbb/\Wbb + 
W bd \ = 1 to within 1%. 




W b + W d = Wbb + Wbd + Wdd + Wdb = Wtot 



(BIO) 
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